#include "undirected_louvain.h"



//logistic mapping
int LogisticMapping(double x0)
{
	double x = x0;
	for (int i = 0; i < 1000; i++) 
	{  // 迭代1000次
		x = 4.0 * x * (1 - x);  // Logistic映射公式
	}
	int random_number = (int)(x * 100);  
	return random_number;
}




void InitGraph::print(void)
{
	cout << "Initial graph G(V: " << m_V << ", E: "
		<< m_E << ") with adjacency list:\n";
	for (int i = 0; i < m_V; i++)
	{
		cout << "(" << i << ") ";
		for (auto v : m_adj[i])
		{
			cout << " -> " << v;
		}
		cout << "\n";
	}
	cout << "\n";
}

vector<vector<int>> InitGraph::getAdj(void)const
{
	return m_adj;
}

int Undirected_Graph::degree(int v)
{
	if (v == 0)
	{
		return degrees[0];
	}
	else
	{
		return degrees[v] - degrees[v - 1];
	}
}

//找到与给定节点 v 相邻的节点和节点v所连接的第一条边的权重地址，
//这样就可以借助该节点的度获取该节点所有连接的边权重。
pair<vector<int>::iterator, vector<double>::iterator>
Undirected_Graph::adjacent(int v)
{
	if (v == 0)
	{
		return make_pair(edges.begin(), weights.begin());   
	}
	else if (weights.size() != 0)
	{
		return make_pair(edges.begin() + degrees[v - 1],       
			weights.begin() + degrees[v - 1]);
	}
	else
	{
		return make_pair(edges.begin() + degrees[v - 1], weights.begin());
	}
}

//判断是否存在自环:在压缩节点时，算法会将模块内部的边权重之和视为自环
double Undirected_Graph::self_loops(int v)
{
	auto pointer = this->adjacent(v);
	for (int i = 0; i < this->degree(v); i++)
	{
		if (*(pointer.first + i) == v)
		{
			if (weights.size() != 0)
				return *(pointer.second + i);  
			else
				return 1.;
		}
	}
	return 0.;  
}

double Undirected_Graph::weighted_degree(int v)
{
	if (weights.size() == 0)
	{
		return this->degree(v);
	}
	else
	{
		auto pointer = this->adjacent(v);
		int degree = this->degree(v);
		double sum = 0.0;
		
		for (int i = 0; i < degree; i++)
		{
			sum += *(pointer.second + i);
		}
		return sum;
	}
}


/*
* 判断节点是否连接，返回两节点之间边的权重
*/
//double Graph::GetNodeConnect(int start, int end)const
//{
//	if (matrix[start][end] != 0)
//	{
//		return matrix[start][end];
//	}
//}


//DEBUG
void Undirected_Graph::print()
{
	cout << "Graph G(V: " << m_V << ", E: " << m_E << ", W: " << m_W << ")\ndegrees: ";
	for (int i = 0; i < m_V; i++)
	{
		cout << degrees[i] << " ";
	}
	cout << "\nedges: ";
	for (int i = 0; i < edges.size(); i++)
	{
		cout << edges[i] << " ";
	}
	cout << "\n";
	cout << "weights: ";
	for (int i = 0; i < weights.size(); i++)
	{
		cout << weights[i] << " ";
	}
	cout << "\n\n";
}




	/*
	* Function:从社区中移除该节点
	* 
	*/
void Undirected_Community::remove(int vertex, int community, double edgesToCommunity)
{
	m_eout[community] -= m_graph.weighted_degree(vertex);
	m_ein[community] -= 2.0 * edgesToCommunity + m_graph.self_loops(vertex);
	m_vertexToCommunityMap[vertex] = -1;
}


void Undirected_Community::insert(int vertex, int community, double edgesToCommunity)
{
	//assert(node >= 0 && node < size);
	m_eout[community] += m_graph.weighted_degree(vertex);
	m_ein[community] += 2.0 * edgesToCommunity + m_graph.self_loops(vertex);
	m_vertexToCommunityMap[vertex] = community;
}


/*
* 模块化公式计算Q
*/
double Undirected_Community::modularity()
{
	double mod = 0.;
	for (int i = 0; i < m_networkSize; i++)
	{
		if (m_eout[i] > 0)
		{
			mod += (m_ein[i] / m_graph.m_W)
				- (m_eout[i] / m_graph.m_W)
				* (m_eout[i] / m_graph.m_W);
		}
	}
	return mod;
}



double Undirected_Community::regularization()
{
	double densitySum = 0.;						
	for (int i = 0; i < m_networkSize; i++)
	{
		if (m_c[i] == 1)
		{
			densitySum++;
		}
		else 
		{
			double den  = 0.5 * m_c[i] * (m_c[i] - 1);
			densitySum += m_ein[i] / den;
		}
	}
	double n = static_cast<double>(m_networkSize);
	double V = static_cast<double>(m_originalV);
	return 0.5 * (static_cast<double>(densitySum) / n - (n / V));	
}


/*
* Function:模块化增益计算（delt_Q）
*/
double Undirected_Community::modularity_gain(int v, 
											 int newCommunity, 
											 double numEdges, 
											 double degree)
{
	double eouts = m_eout[newCommunity];
	double deg   = degree;
	double m     = m_graph.m_W;
	return (numEdges - (eouts * deg) / m); 
}



void Undirected_Community::adj_communities(int vertex)
{
	for (int i = 0; i < m_adjLast; i++)
	{
		m_adjWeights[m_adjPositions[i]] = -1.;
	}

	m_adjLast = 0;
	auto p  = m_graph.adjacent(vertex);
	int deg = m_graph.degree(vertex);

	m_adjPositions[0]				= m_vertexToCommunityMap[vertex];
	m_adjWeights[m_adjPositions[0]] = 0.;
	m_adjLast						= 1;

	for (int i = 0; i < deg; i++) 
	{
		int neighbour		   = *(p.first + i);
		int neighbourCommunity = m_vertexToCommunityMap[neighbour];
		double neighbourWeight = (m_graph.weights.size() == 0) ? 1.0 : *(p.second + i);

		if (neighbour != vertex) 
		{
			if (m_adjWeights[neighbourCommunity] == -1.)
			{
				m_adjWeights[neighbourCommunity] = 0.;
				m_adjPositions[m_adjLast]		 = neighbourCommunity;
				m_adjLast++;
			}
			m_adjWeights[neighbourCommunity] += neighbourWeight;
		}
	}
}

/*
* Function:Louvain算法第一阶段
* 
*/
bool Undirected_Community::phase1(bool dontLimit)
{
	bool hasChanged = false;	
	int moves;
	double newMod	  = modularity();
	double currentMod = newMod; 

	vector<int> randomOrder(m_networkSize);
	for (int i = 0; i < m_networkSize; i++)
	{
		randomOrder[i] = i;  
	}


	for (int i = 0; i < m_networkSize; i++)
	{
		int randPosition = rand() % (m_networkSize - i) + i;
		int tmp			 = randomOrder[i];
		randomOrder[i]	 = randomOrder[randPosition];
		randomOrder[randPosition] = tmp;
	}

	int counter = 0;

	do {
		currentMod = newMod;
		moves = 0;

		for (int i = 0; i < m_networkSize; i++)
		{
			int vertex			= randomOrder[i];
			int vertexCommunity = m_vertexToCommunityMap[vertex];
			double vertexDegree = m_graph.weighted_degree(vertex);

			adj_communities(vertex);
			remove(vertex, vertexCommunity, m_adjWeights[vertexCommunity]);

			int bestCommunity	= vertexCommunity;
			double bestNumEdges = 0.0;
			double bestIncrease = 0.0;		//控制分辨率
			for (int j = 0; j < m_adjLast; j++)
			{
				double increase = modularity_gain(vertex, m_adjPositions[j], m_adjWeights[m_adjPositions[j]], vertexDegree);
				if (increase > bestIncrease)
				{
					bestCommunity = m_adjPositions[j];
					bestNumEdges  = m_adjWeights[m_adjPositions[j]];
					bestIncrease  = increase;
				}
			}
			insert(vertex, bestCommunity, bestNumEdges);

			if (bestCommunity != vertexCommunity)
			{
				moves++;
			}
			counter++;
		}

		newMod = modularity();
		if (moves > 0)
		{
			hasChanged = true;
		}
	} while ((dontLimit || counter < MAX_STEPS) && moves > 0 && newMod > currentMod);

	return hasChanged;
}



Undirected_Graph Undirected_Community::phase2(vector<int>& communitySizes)
{
	vector<int> renumber(m_networkSize, -1);
	for (int v = 0; v < m_networkSize; v++)
	{
		renumber[m_vertexToCommunityMap[v]]++;
	}

	int final = 0;
	for (int i = 0; i < m_networkSize; i++)
	{
		if (renumber[i] != -1) 
		{
			renumber[i] = final;
			final++;
		}
	}

	vector<vector<int>> communityVertices(final);
	for (int v = 0; v < m_networkSize; v++)
	{
		communityVertices[renumber[m_vertexToCommunityMap[v]]].push_back(v);
	}

	Undirected_Graph g;
	int numCommunities = communityVertices.size();
	g.m_V			   = numCommunities;
	communitySizes.resize(numCommunities);
	g.degrees.resize(numCommunities);

	for (int i = 0; i < numCommunities; i++)
	{
		int sum = 0;
		for (int j = 0; j < communityVertices[i].size(); j++) 
		{
			sum += this->m_c[communityVertices[i][j]];
		}
		communitySizes[i] = sum;
	}

	for (int c = 0; c < numCommunities; c++)
	{
		map<int, double> m;
		map<int, double>::iterator it;

		int communitySize = communityVertices[c].size();
		for (int v = 0; v < communitySize; v++)
		{
			auto p  = m_graph.adjacent(communityVertices[c][v]);
			int deg = m_graph.degree(communityVertices[c][v]);
			for (int i = 0; i < deg; i++)
			{
				int neighbour = *(p.first + i);
				int neighbourCommunity = renumber[m_vertexToCommunityMap[neighbour]];
				double neighbourWeight = (m_graph.weights.size() == 0) ? 1 : *(p.second + i);

				it = m.find(neighbourCommunity);
				if (it == m.end())
				{
					m.insert(make_pair(neighbourCommunity, neighbourWeight));
				}
				else
				{
					it->second += neighbourWeight;
				}
			}
		}

		g.degrees[c] = (c == 0) ? m.size() : g.degrees[c - 1] + m.size();
		g.m_E += m.size();

		for (it = m.begin(); it != m.end(); it++)
		{
			g.m_W += it->second;
			g.edges.push_back(it->first);
			g.weights.push_back(it->second);
		}
	}
	return g;
}


void Undirected_Community::partition(vector<int>& result)
{
	vector<int> renumber(m_networkSize, -1);
	for (int v = 0; v < m_networkSize; v++)
	{
		renumber[m_vertexToCommunityMap[v]]++;
	}
	int final = 0;
	for (int i = 0; i < m_networkSize; i++)
	{
		if (renumber[i] != -1)
		{
			renumber[i] = final;
			final++;
		}
	}
	for (int i = 0; i < m_networkSize; i++)
	{
		result[i] = renumber[m_vertexToCommunityMap[i]];
	}
}

//DEBUG
void Undirected_Community::Print()
{
	cout << "Community:\n";
	cout << "Network size = " << m_networkSize << "\n";
	cout << "vertexToCommunityMap: ";
	for (int i = 0; i < m_networkSize; i++)
	{
		cout << m_vertexToCommunityMap[i] << " ";
	}
	cout << "\nein: ";
	for (int i = 0; i < m_networkSize; i++)
	{
		cout << m_ein[i] << " ";
	}
	cout << "\neout: ";
	for (int i = 0; i < m_networkSize; i++)
	{
		cout << m_eout[i] << " ";
	}
	cout << "\nc: ";
	for (int i = 0; i < m_networkSize; i++)
	{
		cout << m_c[i] << " ";
	}
	cout << "\n";
	cout << "modularity: " << modularity() << "\n";
	cout << "reg: " << regularization() << "\n";
}


//获取邻接矩阵或文件中的数据
//pair<unsigned int, unsigned int> GetgraphData( vector<pair<unsigned int, unsigned int>>& edges,
//											   vector<double>& weights, 
//											   Graph_Way way,
//											   const vector<vector<double>>& dsm, 
//											   string filename)
//{
//	pair<unsigned int, unsigned int> edge_node_number;
//	///*******获取图数据***************************//
//	switch (way)
//	{
//	case adjmatrix:
//		edge_node_number = DSMToList(dsm, edges, weights, 1);
//		break;
//	case adjlist://
//		edge_node_number = ReadListData(filename, edges, weights);
//		break;
//	default:
//		break;
//	}
//	return edge_node_number;
//}




map<int, vector<int>> Undirected_Louvain_run(double& q,
											 vector<pair<unsigned int, unsigned int>> edges,
											 vector<double> weights, 
											 pair<unsigned int, unsigned int> dim)
{
	const int V = dim.first; 
	int E = dim.second;            
	InitGraph g(V, E, edges, weights);
	Undirected_Graph graph(g);
	


	double predQ = 0.;

	map<int, vector<int>> best;
	for (int rep = 0; rep < REPS; rep++)
	{		
		Undirected_Graph graph(g);
		vector<int> cs(V, 1);
		Undirected_Community c(graph, V, cs);

		bool hasChanged = true;
		double mod = c.modularity(), newMod;       
		
		double reg = 0.0, newReg;
		vector<float> Qs;
		vector<vector<int>> levels;
		
		int count = 0;
		bool dontLimit = rand() % 3;
		do {			
			// Phase 1: 模块度贪心优化
			hasChanged = c.phase1(dontLimit);
			newMod = c.modularity();
			
			// Phase 2: 社区聚合
			vector<int> result(c.m_networkSize);
			c.partition(result);
			levels.push_back(result);
			vector<int> communitySizes;
			graph = c.phase2(communitySizes);

			c = Undirected_Community(graph, V, communitySizes);
			newReg = c.regularization();
			Qs.push_back(newMod);

			
			mod = newMod;
			reg = newReg;
			
			if (hasChanged)
			{
				count++;
			}
		} while (hasChanged);
		//} while (hasChanged && count < 1);	//利用count获取层次聚类结果

		pair<int, double> maxQ = { 0, 0.0 };
		for (int q = 0; q < Qs.size(); q++)
		{
			if (Qs[q] > maxQ.second)
			{
				maxQ = { q, Qs[q] };
			}
		}

		if (maxQ.second > predQ)
		{
			vector<int> vertexToCommunity(V);
			map<int, vector<int>> communities;

			for (int i = 0; i < V; i++)
			{
				vertexToCommunity[i] = i;
			}


			for (int l = 0; l < levels.size(); l++)
			{
				for (int v = 0; v < V; v++)
				{
					vertexToCommunity[v] = levels[l][vertexToCommunity[v]];
				}					
			}

			for (int v = 0; v < V; v++)
			{
				communities[vertexToCommunity[v]].push_back(v);
			}
			predQ = maxQ.second;
			best = communities;
		}
	} 
	q = predQ;				


	return best;
}


void ListToMatrix(vector<pair<unsigned int, unsigned int>>& edges, 
				  vector<double>& weights, 
				  pair<unsigned int, unsigned int> dim, 
				  vector<vector<double>>& dsm)
{
	int number = 0;
	dsm.resize(dim.first);
	
	//初始化dsm的大小
	for (int i = 0; i < dim.first; i++)
	{
		dsm[i].resize(dim.first);
	}

	for (auto& edg : edges)
	{
		//无向加权图
		if (weights.size() != 0)
		{
			//dsm[edg.first][edg.second] = dsm[edg.second][edg.first] = weights.at(number);
			dsm[edg.first - 1][edg.second - 1] = 
			dsm[edg.second - 1][edg.first - 1] = weights.at(number);
			number++;
		}
		else
		{
			//无向无权图  列表节点标签从1开始，所以这里-1从0开始
			dsm[edg.first - 1][edg.second - 1] = dsm[edg.second - 1][edg.first - 1] = 1.;  
		}
	}
}